Capitulo 4
10
Debemos de utilizar la base de datos Weekly, esta base de datos contiene información del rendimiento porcentual semanal para el índice bursátil S&P 500 entre 1990 y 2010.
## Year Lag1 Lag2 Lag3
## Min. :1990 Min. :-18.1950 Min. :-18.1950 Min. :-18.1950
## 1st Qu.:1995 1st Qu.: -1.1540 1st Qu.: -1.1540 1st Qu.: -1.1580
## Median :2000 Median : 0.2410 Median : 0.2410 Median : 0.2410
## Mean :2000 Mean : 0.1506 Mean : 0.1511 Mean : 0.1472
## 3rd Qu.:2005 3rd Qu.: 1.4050 3rd Qu.: 1.4090 3rd Qu.: 1.4090
## Max. :2010 Max. : 12.0260 Max. : 12.0260 Max. : 12.0260
## Lag4 Lag5 Volume
## Min. :-18.1950 Min. :-18.1950 Min. :0.08747
## 1st Qu.: -1.1580 1st Qu.: -1.1660 1st Qu.:0.33202
## Median : 0.2380 Median : 0.2340 Median :1.00268
## Mean : 0.1458 Mean : 0.1399 Mean :1.57462
## 3rd Qu.: 1.4090 3rd Qu.: 1.4050 3rd Qu.:2.05373
## Max. : 12.0260 Max. : 12.0260 Max. :9.32821
## Today Direction
## Min. :-18.1950 Down:484
## 1st Qu.: -1.1540 Up :605
## Median : 0.2410
## Mean : 0.1499
## 3rd Qu.: 1.4050
## Max. : 12.0260Tenemos una base de datos con 1089 observaciones sobre las siguientes 9 variables.
- Year: El año en que se registró la observación.
- Lag1: Porcentaje de retorno de la semana anterior.
- Lag2: Porcentaje de retorno de 2 semanas anteriores
- Lag3: de retorno de 3 semanas anteriores
- Lag4: Porcentaje de retorno de 4 semanas anteriores
- Lag5: Porcentaje de retorno de 5 semanas anteriores
- Volume: Volumen de acciones negociadas (número promedio de acciones diarias negociadas en miles de millones)
- Today: Porcentaje de retorno para esta semana
- Direction: Un factor con niveles hacia abajo y hacia arriba que indica si el mercado tuvo un rendimiento positivo o negativo en una semana determinada
library(psych)
# Gráfico de dispersión
pairs.panels(Weekly[,-9],
method = "pearson",
hist.col = "#00AFBB",
density = TRUE,
)
Como podemos apreciar, en general no pareciera que hubiesen muchas variables correlacionadas, pero si podemos apreciar una correlación positiva entre el año y la variable Volume. La variable volume hace referencia al volumen de acciones negociadas (número promedio de acciones diarias negociadas en miles de millones), por lo que a medida que pasa el tiempo, podemos apreciar como se aumenta también el volumen de acciones negociadas.
Le ejecutamos a todos los datos, utilizando Direction como la variable respuesta y las variables lag más la variable volume como predictoras, una regresión logística
weeklyrlo <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Weekly, family = "binomial")
summary(weeklyrlo)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = "binomial", data = Weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4Podemos apreciar que la variable lag2, es la única variables significativa.
## Loading required package: lattice## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphapredichos=as.factor(ifelse(test = weeklyrlo$fitted.values > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Weekly$Direction)
matriz<-confusionMatrix(predichos, verdaderos)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Down Up
## Down 54 48
## Up 430 557
##
## Accuracy : 0.5611
## 95% CI : (0.531, 0.5908)
## No Information Rate : 0.5556
## P-Value [Acc > NIR] : 0.369
##
## Kappa : 0.035
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.11157
## Specificity : 0.92066
## Pos Pred Value : 0.52941
## Neg Pred Value : 0.56434
## Prevalence : 0.44444
## Detection Rate : 0.04959
## Detection Prevalence : 0.09366
## Balanced Accuracy : 0.51612
##
## 'Positive' Class : Down
## | Reference | ||
|---|---|---|
| Predicted | Down | Up |
| Down | A | B |
| Up | C | D |
- Lo primero que podemos apreciar es que la precisión del modelo es de un 56.11% (Accuracy).
- También podemos ver que la sensitividad (\(Sensitivity = \frac{A}{A + C}\)), es de solamente un 11.15%, lo cual quiere decir que en semanas en donde el mercado es bajista.
- La especificidad (\(Specificity = \frac{D}{B + D}\)), es de un 92.06%, es la precisión a la hora de predecir semanas alcistas
weeklyrlo2 <- glm(Direction ~ Lag2, data = Weekly, family = "binomial", subset = w_entrenamiento)
summary(weeklyrlo2)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = Weekly,
## subset = w_entrenamiento)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4Weekly.20092010 <- Weekly[!w_entrenamiento, ]
Direction.20092010 <- Direction[!w_entrenamiento]
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Down Up
## Down 9 5
## Up 34 56
##
## Accuracy : 0.625
## 95% CI : (0.5247, 0.718)
## No Information Rate : 0.5865
## P-Value [Acc > NIR] : 0.2439
##
## Kappa : 0.1414
##
## Mcnemar's Test P-Value : 7.34e-06
##
## Sensitivity : 0.20930
## Specificity : 0.91803
## Pos Pred Value : 0.64286
## Neg Pred Value : 0.62222
## Prevalence : 0.41346
## Detection Rate : 0.08654
## Detection Prevalence : 0.13462
## Balanced Accuracy : 0.56367
##
## 'Positive' Class : Down
## - Podemos apreciar que la precisión del modelo es de un 62.50%
- Cuando la semana del mercado es bajista, la precisión es de un 20.93%
- Cuando la semana del mercado es alcista, la precisión es de un 91.80%
## Call:
## lda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162pred.lda <- predict(ldapre, Weekly.20092010)
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Down Up
## Down 9 5
## Up 34 56
##
## Accuracy : 0.625
## 95% CI : (0.5247, 0.718)
## No Information Rate : 0.5865
## P-Value [Acc > NIR] : 0.2439
##
## Kappa : 0.1414
##
## Mcnemar's Test P-Value : 7.34e-06
##
## Sensitivity : 0.20930
## Specificity : 0.91803
## Pos Pred Value : 0.64286
## Neg Pred Value : 0.62222
## Prevalence : 0.41346
## Detection Rate : 0.08654
## Detection Prevalence : 0.13462
## Balanced Accuracy : 0.56367
##
## 'Positive' Class : Down
## - La precisión de este modelo es de 62.50%
- Cuando en la semana el mercado es bajista la precisión es de 20.93%
- Cuando en la semana el mercado es alcista la precisión es de 91.80%
## Call:
## qda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Down Up
## Down 0 0
## Up 43 61
##
## Accuracy : 0.5865
## 95% CI : (0.4858, 0.6823)
## No Information Rate : 0.5865
## P-Value [Acc > NIR] : 0.5419
##
## Kappa : 0
##
## Mcnemar's Test P-Value : 1.504e-10
##
## Sensitivity : 0.0000
## Specificity : 1.0000
## Pos Pred Value : NaN
## Neg Pred Value : 0.5865
## Prevalence : 0.4135
## Detection Rate : 0.0000
## Detection Prevalence : 0.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : Down
## - En este caso la predicción del modelo es de un 58.65%
- La precisión cuando es mercado es bajista es de un 0%
- La precisión cuando el mercado es alcista es de un 100%
library(class)
entrenamiento.x <- as.matrix(Lag2[w_entrenamiento])
prueba.x <- as.matrix(Lag2[!w_entrenamiento])
entrenamiento.Direction <- Direction[w_entrenamiento]
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 1)
confusionMatrix(pred.knn, Direction.20092010)## Confusion Matrix and Statistics
##
## Reference
## Prediction Down Up
## Down 21 30
## Up 22 31
##
## Accuracy : 0.5
## 95% CI : (0.4003, 0.5997)
## No Information Rate : 0.5865
## P-Value [Acc > NIR] : 0.9700
##
## Kappa : -0.0033
##
## Mcnemar's Test P-Value : 0.3317
##
## Sensitivity : 0.4884
## Specificity : 0.5082
## Pos Pred Value : 0.4118
## Neg Pred Value : 0.5849
## Prevalence : 0.4135
## Detection Rate : 0.2019
## Detection Prevalence : 0.4904
## Balanced Accuracy : 0.4983
##
## 'Positive' Class : Down
## - La precisión del modelo es de un 50%
- Cuando el mercado es bajista la precisión es de un 48.84%
- Cuando el mercado es alcista la precisión es de un 50.82%
- Si analizamos las precisiones de los modelos obtuvimos los siguientes resultados:
| GLN | LDA | QDA | KNN | |
|---|---|---|---|---|
| Precisión | 62.5 | 62.5 | 58.65 | 50 |
| Sensitividad | 20.93 | 20.93 | 0 | 48.84 |
| Especificidad | 91.80 | 91.80 | 100 | 50.82 |
Por lo que podemos concluir que tanto el modelo GLN como LDA fueron los que mejor resultados dieron.
A continuación evaluamos KNN con los valores de K = 10, 20, 50 y 100
# Para K = 10
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 10)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table## Reference
## Prediction Down Up
## Down 16 20
## Up 27 41## [1] "Accuracy: 0.5481 Sensitivity: 0.3721 Specificity: 0.6721"# Para k = 20
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 20)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table## Reference
## Prediction Down Up
## Down 23 21
## Up 20 40## [1] "Accuracy: 0.6058 Sensitivity: 0.5349 Specificity: 0.6557"# Para k = 50
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 50)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table## Reference
## Prediction Down Up
## Down 20 22
## Up 23 39## [1] "Accuracy: 0.5673 Sensitivity: 0.4651 Specificity: 0.6393"# Para k = 50
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 100)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table## Reference
## Prediction Down Up
## Down 10 11
## Up 33 50## [1] "Accuracy: 0.5769 Sensitivity: 0.2326 Specificity: 0.8197"Ahora probaremos una regresión logística con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume
# para volume
weeklyrlo2 <- glm(Direction ~ Volume, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 31 47
## Up 12 14## [1] "Accuracy: 0.4327 Sensitivity: 0.7209 Specificity: 0.2295"# para lag2 + lag3
weeklyrlo2 <- glm(Direction ~ Lag2 + Lag3, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 8 4
## Up 35 57## [1] "Accuracy: 0.625 Sensitivity: 0.186 Specificity: 0.9344"# para lag1 + lag2
weeklyrlo2 <- glm(Direction ~ Lag1 + Lag2, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 7 8
## Up 36 53## [1] "Accuracy: 0.5769 Sensitivity: 0.1628 Specificity: 0.8689"# para lag1 + lag2 + volume
weeklyrlo2 <- glm(Direction ~ Lag1 + Lag2 + Volume, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 27 33
## Up 16 28## [1] "Accuracy: 0.5288 Sensitivity: 0.6279 Specificity: 0.459"Ahora probaremos una LDA con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume
# Para volume
ldapre <- lda(Direction ~ Volume, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 33 47
## Up 10 14## [1] "Accuracy: 0.4519 Sensitivity: 0.7674 Specificity: 0.2295"# Para lag2 + lag3
ldapre <- lda(Direction ~ Lag2 + Lag3, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 8 4
## Up 35 57## [1] "Accuracy: 0.625 Sensitivity: 0.186 Specificity: 0.9344"# Para lag1 + lag2
ldapre <- lda(Direction ~ Lag1 + Lag2, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 7 8
## Up 36 53## [1] "Accuracy: 0.5769 Sensitivity: 0.1628 Specificity: 0.8689"# Para lag1 + lag2 + volume
ldapre <- lda(Direction ~ Lag1 + Lag2 + Volume, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 27 33
## Up 16 28## [1] "Accuracy: 0.5288 Sensitivity: 0.6279 Specificity: 0.459"Ahora probaremos una QDA con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume
# para volume
qdapre <- qda(Direction ~ Volume, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 43 59
## Up 0 2## [1] "Accuracy: 0.4327 Sensitivity: 1 Specificity: 0.0328"# para lag2 + lag3
qdapre <- qda(Direction ~ Lag2 + Lag3, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 4 2
## Up 39 59## [1] "Accuracy: 0.6058 Sensitivity: 0.093 Specificity: 0.9672"# para lag1 + lag2
qdapre <- qda(Direction ~ Lag1 + Lag2, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 7 10
## Up 36 51## [1] "Accuracy: 0.5577 Sensitivity: 0.1628 Specificity: 0.8361"# para lag1 + lag2 + volume
qdapre <- qda(Direction ~ Lag1 + Lag2 + Volume, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table## Reference
## Prediction Down Up
## Down 31 44
## Up 12 17## [1] "Accuracy: 0.4615 Sensitivity: 0.7209 Specificity: 0.2787"11
Auto <- Auto[c(-10,-11,-12,-13,-14,-15,-16,-17,-18)]
mpg01 <- rep(0, length(Auto$mpg))
mpg01[Auto$mpg > median(Auto$mpg)] <- 1
Auto <- data.frame(Auto, mpg01)## 'data.frame': 392 obs. of 10 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
## $ mpg01 : num 0 0 0 0 0 0 0 0 0 0 ...# Sabiendo que name es un factor, lo excluimos del análisis
pairs.panels(Auto[,-9],
method = "pearson",
hist.col = "#00AFBB",
density = TRUE,
)
Podemos apreciar que existe una correlación negativa con las variables cylinders, displacement, horsepower y weight
## Call:
## lda(mpg01 ~ displacement + weight + cylinders + horsepower, data = train)
##
## Prior probabilities of groups:
## 0 1
## 0.5 0.5
##
## Group means:
## displacement weight cylinders horsepower
## 0 271.4830 3601.537 6.789116 129.09524
## 1 118.0782 2352.027 4.217687 79.47619
##
## Coefficients of linear discriminants:
## LD1
## displacement -9.297213e-04
## weight -1.046022e-03
## cylinders -3.826930e-01
## horsepower 7.100215e-05pred.lda <- predict(autos.lda, test, type="response")
matriz <- confusionMatrix(pred.lda$class, as.factor(test$mpg01))
matriz$table## Reference
## Prediction 0 1
## 0 42 1
## 1 7 48## [1] "Accuracy: 0.9184 Sensitivity: 0.8571 Specificity: 0.9796"Podemos apreciar que la tasa de error es de un 10.2%
## Call:
## qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = train)
##
## Prior probabilities of groups:
## 0 1
## 0.5 0.5
##
## Group means:
## cylinders weight displacement horsepower
## 0 6.789116 3601.537 271.4830 129.09524
## 1 4.217687 2352.027 118.0782 79.47619pred.qda <- predict(autos.qda, test)
matriz <- confusionMatrix(pred.qda$class, as.factor(test$mpg01))
matriz$table## Reference
## Prediction 0 1
## 0 43 3
## 1 6 46## [1] "Accuracy: 0.9082 Sensitivity: 0.8776 Specificity: 0.9388"Como podemos ver la tasa de error es del 11.22%
autos.glm <- glm(mpg01 ~ cylinders + weight + displacement + horsepower, data = train, family = binomial)
summary(autos.glm)##
## Call:
## glm(formula = mpg01 ~ cylinders + weight + displacement + horsepower,
## family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3670 -0.1738 0.0466 0.3810 3.1933
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.7659288 1.9131492 6.150 7.75e-10 ***
## cylinders -0.0161152 0.3729110 -0.043 0.96553
## weight -0.0020281 0.0008052 -2.519 0.01177 *
## displacement -0.0092265 0.0090618 -1.018 0.30859
## horsepower -0.0447260 0.0155085 -2.884 0.00393 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 407.57 on 293 degrees of freedom
## Residual deviance: 163.52 on 289 degrees of freedom
## AIC: 173.52
##
## Number of Fisher Scoring iterations: 7p <- predict(autos.glm, test, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = 1, no = 0))
verdaderos=as.factor(test$mpg01)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 43 1
## 1 6 48## [1] "Accuracy: 0.9286 Sensitivity: 0.8776 Specificity: 0.9796"La tasa de error es de 11.22%
train.X2 = cbind(train$displacement, train$weight, train$cylinders, train$year)
test.X2 = cbind(test$displacement, test$weight, test$cylinders, test$year)
train.Y2 = cbind(train$mpg01)# Para K = 10
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 10)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table## Reference
## Prediction 0 1
## 0 44 2
## 1 5 47## [1] "Accuracy: 0.9286 Sensitivity: 0.898 Specificity: 0.9592"La tasa de error el de 10.2%
# Para K = 50
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 50)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table## Reference
## Prediction 0 1
## 0 43 3
## 1 6 46## [1] "Accuracy: 0.9082 Sensitivity: 0.8776 Specificity: 0.9388"La tasa de error es del 10.2%
# Para K = 100
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 100)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table## Reference
## Prediction 0 1
## 0 42 3
## 1 7 46## [1] "Accuracy: 0.898 Sensitivity: 0.8571 Specificity: 0.9388"La tasa de error es del 10.2%
Podemos ver que a partir de k = 10, la tasa de error se mantiene.
12
## [1] 8## [1] 6561## [1] 1000## [1] 2.2518e+15## [1] 2248091x <- 1:10
plot(x, Power3(x, 2), log = "xy", xlab = "Log de x", ylab = "Log de x^2", main = "Log de x^2 vs Log de x")
13
Se requiere usar la base de datos “Boston”, la cual a continuación podemos apreciar las variables que posee:
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...Ahora calcularemos la mediana de la variiable crim para crear una nueva variable, la cual será binaria en donde 0 = crim <= mediana ó 1 = crim > mediana
mediana_crim <- rep(0, length(Boston$crim))
mediana_crim[Boston$crim > median(Boston$crim)] <- 1
Boston$mediana_crim = mediana_crim
head(Boston)## crim zn indus chas nox rm age dis rad tax ptratio black
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12
## lstat medv mediana_crim
## 1 4.98 24.0 0
## 2 9.14 21.6 0
## 3 4.03 34.7 0
## 4 2.94 33.4 0
## 5 5.33 36.2 0
## 6 5.21 28.7 0Ahora construiremos nuestros datos de entrenamiento y de prueba
set.seed(123)
train_index <- sample(1:nrow(Boston), 0.8 * nrow(Boston))
test_index <- setdiff(1:nrow(Boston), train_index)Ahora para entender un poco mejor las variables realizaremos un gráfico de dispersión y de correlaciones
Podemos apreciar que respecto a la variable mediana_crim, las variables nox, rad, dis, age, tax e indus son las que más están correlacionadas. Para este ejercicio se utilizarán modelos partiendo de una variable (la de mayor correlación) y se irán agregando nuevas variables para determinar cual es el mejor modelo para cada ténica.
Comenzando con modelos de regresión logística:
- Modelo con nox como única variable predictora
glm_13.fit <- glm(mediana_crim ~ nox, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
verdaderos=as.factor(test[,15])
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 40 5
## 1 9 48## [1] "Accuracy: 0.8627 Sensitivity: 0.8163 Specificity: 0.9057"- Modelo con nox y rad como únicas variables predictoras
glm_13.fit <- glm(mediana_crim ~ nox + rad, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 46 7
## 1 3 46## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"- Modelo con nox, rad y dis como únicas variables predictoras
glm_13.fit <- glm(mediana_crim ~ nox + rad + dis, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 46 7
## 1 3 46## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"Como podemos apreciar en los anteriores modelos, con las variables nox y rad como únicas predictoras se obtuvo el mejor modelo, cuando agregamos más variables predictoras no se mejoran los resultados: Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679
Ahora seguimos con los modelos LDA
- Modelo con nox como única variable predictora
ldapre <- lda(mediana_crim ~ nox, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 42 11
## 1 7 42## [1] "Accuracy: 0.8235 Sensitivity: 0.8571 Specificity: 0.7925"- Modelo con nox y rad como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 47 12
## 1 2 41## [1] "Accuracy: 0.8627 Sensitivity: 0.9592 Specificity: 0.7736"- Modelo con nox, rad y dis como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad + dis, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 45 11
## 1 4 42## [1] "Accuracy: 0.8529 Sensitivity: 0.9184 Specificity: 0.7925"- Modelo con nox, rad, dis y age como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad + dis + age, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table## Reference
## Prediction 0 1
## 0 45 13
## 1 4 40## [1] "Accuracy: 0.8333 Sensitivity: 0.9184 Specificity: 0.7547"Podemos ver que aumentando el número de variables no presenta un mejora en la precisión de nuestro modelo, de lo que se probaron el mejor fue el que se usó con las variables nox y rad con los siguientes parámetros: Accuracy: 0.8627 Sensitivity: 0.9592 Specificity: 0.7736
Modelos KNN
modelos_knn <- function(x_train, x_test, y_train, y_test, nk){
pred.knn <- knn(x_train, x_test, y_train, k = nk)
z <-confusionMatrix(pred.knn, y_test)
print(paste0("K = ", nk))
imprimirASS(z)
}entrenamiento.x <- as.matrix(train[,-15])
prueba.x <- as.matrix(test[,-15])
entrenamiento.y <- as.factor(train[,15])
prueba.y <- as.factor(test[,15])## [1] "K = 2"
## [1] "Accuracy: 0.9118 Sensitivity: 0.9184 Specificity: 0.9057"
## [1] "K = 3"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9184 Specificity: 0.9245"
## [1] "K = 4"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9388 Specificity: 0.9245"
## [1] "K = 5"
## [1] "Accuracy: 0.951 Sensitivity: 0.9796 Specificity: 0.9245"
## [1] "K = 6"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 7"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 8"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 9"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 10"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 11"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 12"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9592 Specificity: 0.8868"
## [1] "K = 13"
## [1] "Accuracy: 0.9118 Sensitivity: 0.9388 Specificity: 0.8868"
## [1] "K = 14"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9388 Specificity: 0.9057"
## [1] "K = 15"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 16"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 17"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 18"
## [1] "Accuracy: 0.8922 Sensitivity: 0.9184 Specificity: 0.8679"
## [1] "K = 19"
## [1] "Accuracy: 0.8725 Sensitivity: 0.9184 Specificity: 0.8302"
## [1] "K = 20"
## [1] "Accuracy: 0.8824 Sensitivity: 0.9184 Specificity: 0.8491"Podemos ver con un K = 5 se obtuvo la mayor precisión, lo cual también determina, que éste es el mejor modelo de todos los que se probaron en el ejercicio
Capitulo 8
7
## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margin## The following object is masked from 'package:psych':
##
## outlierset.seed(1)
train <- sample(1:nrow(Boston), 400)
Boston.train <- Boston[train, -14]
Boston.test <- Boston[-train, -14]
Y.train <- Boston[train, 14]
Y.test <- Boston[-train, 14]
boston1 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = ncol(Boston) - 1, ntree = 500)
boston2 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = (ncol(Boston) - 1) / 2, ntree = 500)
boston3 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = sqrt(ncol(Boston) - 1), ntree = 500)
valores <- data.frame(c(1:500))
names(valores) <- "arboles"
valores["MSE1"] <- boston1$test$mse
valores["MSE2"] <- boston2$test$mse
valores["MSE3"] <- boston3$test$mse
library(ggplot2)
library(reshape2)
dd = melt(valores, id=c("arboles"))
theme_set(theme_bw())
ggplot(dd) + geom_line(aes(x=arboles, y=value, color=variable)) +
scale_color_manual(name = "m = ",
labels = c("p","p/2","\u221Ap"),
values=c("green","red","blue")) +
labs(x = "Número de árboles",
y = "Error de clasificación") +
theme(legend.position="top")
Como podemos apreciar en el gráfico anterior con solamente un árbol el erro tiende a ser muy alto, pero a medida que se aumenta el número de árboles, se estabiliza este error, aproximadamente a partir de 100 árboles, este comportamiento se presenta en general con los tres valores de m. Respecto a con cual m se presenta el menor error, podemos ver que en mi caso fue con \(m=p\), aunque si cambiamos la semilla este resultado puede cambiar.
8
train <- sample(1:nrow(Carseats), 300)
Carseats.train <- Carseats[train, ]
Carseats.test <- Carseats[-train, ]library(rpart)
library(rpart.plot)
tree <- rpart(Sales~., data=Carseats.train)
rpart.plot(tree, box.palette="RdBu", nn=TRUE)
## [1] 4.197228Como podemos apreciar el error MSE de la prueba es de 4.197228
##
## Regression tree:
## rpart(formula = Sales ~ ., data = Carseats.train)
##
## Variables actually used in tree construction:
## [1] Advertising Age CompPrice Income Population Price
## [7] ShelveLoc
##
## Root node error: 2247.1/300 = 7.4904
##
## n= 300
##
## CP nsplit rel error xerror xstd
## 1 0.211723 0 1.00000 1.01804 0.082943
## 2 0.123518 1 0.78828 0.80759 0.065671
## 3 0.057792 2 0.66476 0.71205 0.057322
## 4 0.039590 3 0.60697 0.68368 0.053027
## 5 0.032717 4 0.56738 0.67813 0.053775
## 6 0.026549 5 0.53466 0.66341 0.050763
## 7 0.021459 6 0.50811 0.62629 0.048310
## 8 0.018375 7 0.48665 0.63310 0.050340
## 9 0.017637 9 0.44990 0.64361 0.051830
## 10 0.017500 10 0.43226 0.64903 0.051778
## 11 0.015287 11 0.41476 0.62861 0.050278
## 12 0.014037 13 0.38419 0.62867 0.050364
## 13 0.013272 14 0.37015 0.64636 0.051480
## 14 0.012158 15 0.35688 0.64895 0.051655
## 15 0.011647 16 0.34472 0.64894 0.051919
## 16 0.010000 17 0.33308 0.64854 0.051978
## 7
## 7## [1] 0.02145892Con la validación cruzada obtenemos un tamaño de arbol igual a 7 y un cp de 0.02145892
## [1] 4.54531Podemos ver que se aumenta el MSE en nuestro caso después de la poda
bag.carseats <- randomForest(Sales ~ ., data = Carseats.train, mtry = (ncol(Carseats) - 1), importance = TRUE)
yhat.bag <- predict(bag.carseats, newdata = Carseats.test)
mean((yhat.bag - Carseats.test$Sales)^2)## [1] 2.381242Podemos ver como se mejoró a 2.38 el MSE con este enfoque.
## %IncMSE IncNodePurity
## CompPrice 38.928557 272.987960
## Income 11.725449 122.545479
## Advertising 14.937292 128.773274
## Population 2.518008 95.322579
## Price 71.644649 716.074862
## ShelveLoc 69.101715 596.025736
## Age 17.537409 167.453663
## Education 2.793162 64.343849
## Urban -2.073262 8.985155
## US 4.206136 9.651551Y con esto podemos apreciar que Price, ShelveLoc y CompPrice son las variables predictoras más importantes en este análisis
bag.carseats <- randomForest(Sales ~ ., data = Carseats.train, mtry = sqrt((ncol(Carseats) - 1)), importance = TRUE)
yhat.bag <- predict(bag.carseats, newdata = Carseats.test)
mean((yhat.bag - Carseats.test$Sales)^2)## [1] 2.80204En este caso tenemos un MSE de 2.8
## %IncMSE IncNodePurity
## CompPrice 16.2755940 216.47915
## Income 7.8217396 184.80927
## Advertising 14.1807157 183.03683
## Population 0.6892522 155.21205
## Price 42.0256227 551.81898
## ShelveLoc 47.8554474 473.99364
## Age 12.1305637 213.32128
## Education 2.8662179 91.35868
## Urban -2.3838218 18.20191
## US 3.4288791 26.69131También podemos ver que las variables más importantes siguen siendo las mismas que las que se declararon en el punto anterior.
9
##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "SpecialCH" "PriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.707 = 559.3 / 791
## Misclassification error rate: 0.1512 = 121 / 800Como podemos apreciar, el arbol tiene 9 nodos terminales y un error de entrenamiento de 0.1512. También podemos ver que las variables usadas en la construcción del árbol son LoyalCH, SalePriceMM, specialCH y PriceDiff
## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1065.00 CH ( 0.61625 0.38375 )
## 2) LoyalCH < 0.5036 351 424.90 MM ( 0.29345 0.70655 )
## 4) LoyalCH < 0.264232 151 102.10 MM ( 0.10596 0.89404 )
## 8) LoyalCH < 0.051325 57 0.00 MM ( 0.00000 1.00000 ) *
## 9) LoyalCH > 0.051325 94 85.77 MM ( 0.17021 0.82979 ) *
## 5) LoyalCH > 0.264232 200 273.90 MM ( 0.43500 0.56500 )
## 10) SalePriceMM < 2.04 109 126.30 MM ( 0.26606 0.73394 )
## 20) SpecialCH < 0.5 83 78.43 MM ( 0.18072 0.81928 ) *
## 21) SpecialCH > 0.5 26 35.89 CH ( 0.53846 0.46154 ) *
## 11) SalePriceMM > 2.04 91 119.20 CH ( 0.63736 0.36264 ) *
## 3) LoyalCH > 0.5036 449 349.40 CH ( 0.86860 0.13140 )
## 6) LoyalCH < 0.764572 188 217.80 CH ( 0.73404 0.26596 )
## 12) PriceDiff < 0.265 113 153.40 CH ( 0.58407 0.41593 )
## 24) PriceDiff < -0.165 34 41.19 MM ( 0.29412 0.70588 ) *
## 25) PriceDiff > -0.165 79 95.30 CH ( 0.70886 0.29114 ) *
## 13) PriceDiff > 0.265 75 25.19 CH ( 0.96000 0.04000 ) *
## 7) LoyalCH > 0.764572 261 78.30 CH ( 0.96552 0.03448 ) *Se escoge el nodo terminal 7, en este podemos apreciar que su creterio de división es LoyalCH > 0.764572, éste tiene 261 observaciones en este nodo con una desviación de 78.30. La predicción en este nodo es Sales = CH, también podemos ver que el 96.552% de las observaciones toman el valor de CH, mientras que el 3.448% toman el valor de MM.
En este caso podemos apreciar que la variable principal es LoyalCH, los nodos que nacen de este también usan LoyalCH para dividir el arbol. Podemos ver que en el nodo principal si LoyalCH es mayor o igual a 0.5036, sólo se usa PriceDiff como otra variable, mientras en caso contrario se utilizan las variables salePriceMM y SpecialCH.
tree.pred <- predict(tree.oj, OJ.test, type = "class")
matriz<-confusionMatrix(tree.pred, OJ.test$Purchase)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 144 34
## MM 16 76
##
## Accuracy : 0.8148
## 95% CI : (0.7633, 0.8593)
## No Information Rate : 0.5926
## P-Value [Acc > NIR] : 4.577e-15
##
## Kappa : 0.6064
##
## Mcnemar's Test P-Value : 0.01621
##
## Sensitivity : 0.9000
## Specificity : 0.6909
## Pos Pred Value : 0.8090
## Neg Pred Value : 0.8261
## Prevalence : 0.5926
## Detection Rate : 0.5333
## Detection Prevalence : 0.6593
## Balanced Accuracy : 0.7955
##
## 'Positive' Class : CH
## Podemos ver que la precisión del modelo es de 81.48%, siendo mucho más preciso para predecir CH con un 90% mientras que para predecir MM sólo es de un 69.09%.
## $size
## [1] 9 8 7 4 2 1
##
## $dev
## [1] 135 135 140 151 172 303
##
## $k
## [1] -Inf 0.000000 2.000000 4.666667 12.500000 145.000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"valores <- data.frame(cv.oj$size, cv.oj$dev)
names(valores) <- c("Size", "Dev")
library(ggplot2)
library(reshape2)
theme_set(theme_bw())
ggplot(data = valores, aes(x=Size, y=Dev)) +
geom_line(color="red", linetype = "dashed") +
geom_point() +
scale_x_discrete(limits=c(1:9))
El tamaño de árbol con menor error es de 7
##
## Classification tree:
## snip.tree(tree = tree.oj, nodes = c(4L, 10L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "PriceDiff"
## Number of terminal nodes: 7
## Residual mean deviance: 0.7409 = 587.5 / 793
## Misclassification error rate: 0.1538 = 123 / 800##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "SpecialCH" "PriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.707 = 559.3 / 791
## Misclassification error rate: 0.1512 = 121 / 800Vemos que practicamente no existe diferencia entre los dos árboles, siendo levemenete mayor en el arbol podado
prune.pred <- predict(prune.oj, OJ.test, type = "class")
matriz<-confusionMatrix(prune.pred, OJ.test$Purchase)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 141 34
## MM 19 76
##
## Accuracy : 0.8037
## 95% CI : (0.7512, 0.8494)
## No Information Rate : 0.5926
## P-Value [Acc > NIR] : 1.152e-13
##
## Kappa : 0.5846
##
## Mcnemar's Test P-Value : 0.05447
##
## Sensitivity : 0.8812
## Specificity : 0.6909
## Pos Pred Value : 0.8057
## Neg Pred Value : 0.8000
## Prevalence : 0.5926
## Detection Rate : 0.5222
## Detection Prevalence : 0.6481
## Balanced Accuracy : 0.7861
##
## 'Positive' Class : CH
## Vemos que bajó la precisión del modelo con respecto al árbol sin podar.
10
## Loaded gbm 2.1.5lambdas = seq(0.001, 0.3, 0.005)
mse <- rep(0, length(lambdas))
for (i in 1:length(lambdas)) {
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
pred.train <- predict(boost.hitters, Hitters.train, n.trees = 1000)
mse[i] <- mean((pred.train - Hitters.train$Salary)^2)
}library(ggplot2)
mse_graph <- data.frame(lambdas)
mse_graph["MSE"]<- mse
names(mse_graph) <- c("Lambdas", "MSE")
theme_set(theme_bw())
ggplot(data = mse_graph, aes(x = Lambdas, y = MSE)) +
geom_line(color="red", linetype = "dashed") +
geom_point() 
mse2 <- rep(NA, length(lambdas))
for (i in 1:length(lambdas)) {
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
yhat <- predict(boost.hitters, Hitters.test, n.trees = 1000)
mse2[i] <- mean((yhat - Hitters.test$Salary)^2)
}library(ggplot2)
mse_graph2 <- data.frame(lambdas)
mse_graph2["MSE"]<- mse2
names(mse_graph2) <- c("Lambdas", "MSE")
theme_set(theme_bw())
ggplot(data = mse_graph2, aes(x = Lambdas, y = MSE)) +
geom_line(color="red", linetype = "dashed") +
geom_point() 
## [1] 0.2475835## [1] 0.166lm810 <- lm(Salary ~ ., data = Hitters.train)
pred810 <- predict(lm810, Hitters.test)
mean((pred810 - Hitters.test$Salary)^2)## [1] 0.4917959##
## Attaching package: 'pls'## The following object is masked from 'package:caret':
##
## R2## The following object is masked from 'package:stats':
##
## loadingspcr.fit=pcr(Salary~., data=Hitters.train , scale=TRUE , validation ="CV")
validationplot(pcr.fit ,val.type="MSEP")
## [1] 0.4661183Podemos ver que el MSE por boosting es menor que el dado por la regresión lineal y componentes principales.
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[which.min(mse2)])
summary(boost.hitters)
## var rel.inf
## CAtBat CAtBat 17.7715143
## CWalks CWalks 9.0983419
## CRuns CRuns 8.1091446
## PutOuts PutOuts 7.8019290
## Walks Walks 7.3620230
## CHmRun CHmRun 6.4505765
## Assists Assists 5.6779659
## Hits Hits 5.0562630
## RBI RBI 4.9022635
## Years Years 4.8919157
## AtBat AtBat 4.3145631
## CHits CHits 4.0969282
## CRBI CRBI 3.8001772
## Runs Runs 3.6363269
## HmRun HmRun 3.3847743
## Errors Errors 2.3407756
## NewLeague NewLeague 0.5117111
## Division Division 0.4565112
## League League 0.3362948La variable CatBat es la más importante.
bag.hitters <- randomForest(Salary ~ ., data = Hitters.train, mtry = 19)
yhat.bag <- predict(bag.hitters, newdata = Hitters.test)
mean((yhat.bag - Hitters.test$Salary)^2)## [1] 0.22971Podemos ver que el MSE fue de 0.007 lo cual significa que es mucho mejor que los probados hasta ahora.
11
train <- 1:1000
Caravan$Purchase <- ifelse(Caravan$Purchase == "Yes", 1, 0)
Caravan.train <- Caravan[train, ]
Caravan.test <- Caravan[-train, ]bc <- gbm(Purchase ~ ., data = Caravan.train, distribution = "gaussian", n.trees = 1000, shrinkage = 0.01)## Warning in gbm.fit(x = x, y = y, offset = offset, distribution =
## distribution, : variable 50: PVRAAUT has no variation.## Warning in gbm.fit(x = x, y = y, offset = offset, distribution =
## distribution, : variable 71: AVRAAUT has no variation.
## var rel.inf
## PPERSAUT PPERSAUT 14.6627387
## MKOOPKLA MKOOPKLA 11.5135457
## MOPLHOOG MOPLHOOG 8.1166824
## MBERMIDD MBERMIDD 5.5659055
## PBRAND PBRAND 4.8848882
## MGODGE MGODGE 4.3094073
## ABRAND ABRAND 4.0899213
## MINK3045 MINK3045 3.7799757
## MAUT1 MAUT1 3.0595852
## MOSTYPE MOSTYPE 2.6279165
## PWAPART PWAPART 2.6176353
## MBERARBG MBERARBG 2.1522326
## MGODPR MGODPR 2.1053373
## MAUT2 MAUT2 2.0837180
## MSKC MSKC 2.0548149
## MSKA MSKA 1.8052960
## MSKB1 MSKB1 1.7245921
## PBYSTAND PBYSTAND 1.5426299
## MBERHOOG MBERHOOG 1.4526150
## MFWEKIND MFWEKIND 1.4470476
## MINK7512 MINK7512 1.3601434
## MGODRK MGODRK 1.3247916
## MGODOV MGODOV 1.2359721
## MOPLMIDD MOPLMIDD 1.1189140
## MINK4575 MINK4575 1.0021078
## MRELGE MRELGE 0.9752463
## MRELOV MRELOV 0.9168410
## MAUT0 MAUT0 0.8624586
## MINKGEM MINKGEM 0.8250010
## MFGEKIND MFGEKIND 0.7750144
## APERSAUT APERSAUT 0.7452557
## MBERBOER MBERBOER 0.6269871
## MINKM30 MINKM30 0.6230703
## MOSHOOFD MOSHOOFD 0.6128366
## MZFONDS MZFONDS 0.5781958
## PMOTSCO PMOTSCO 0.5172908
## MHHUUR MHHUUR 0.4691845
## MHKOOP MHKOOP 0.4467587
## MZPART MZPART 0.4351830
## MBERARBO MBERARBO 0.4270172
## MGEMLEEF MGEMLEEF 0.3961990
## MGEMOMV MGEMOMV 0.3886987
## MSKD MSKD 0.3409863
## PLEVEN PLEVEN 0.3215526
## MSKB2 MSKB2 0.2531789
## MINK123M MINK123M 0.2416957
## MOPLLAAG MOPLLAAG 0.1897649
## MRELSA MRELSA 0.1837443
## MBERZELF MBERZELF 0.1086964
## MFALLEEN MFALLEEN 0.1007280
## MAANTHUI MAANTHUI 0.0000000
## PWABEDR PWABEDR 0.0000000
## PWALAND PWALAND 0.0000000
## PBESAUT PBESAUT 0.0000000
## PVRAAUT PVRAAUT 0.0000000
## PAANHANG PAANHANG 0.0000000
## PTRACTOR PTRACTOR 0.0000000
## PWERKT PWERKT 0.0000000
## PBROM PBROM 0.0000000
## PPERSONG PPERSONG 0.0000000
## PGEZONG PGEZONG 0.0000000
## PWAOREG PWAOREG 0.0000000
## PZEILPL PZEILPL 0.0000000
## PPLEZIER PPLEZIER 0.0000000
## PFIETS PFIETS 0.0000000
## PINBOED PINBOED 0.0000000
## AWAPART AWAPART 0.0000000
## AWABEDR AWABEDR 0.0000000
## AWALAND AWALAND 0.0000000
## ABESAUT ABESAUT 0.0000000
## AMOTSCO AMOTSCO 0.0000000
## AVRAAUT AVRAAUT 0.0000000
## AAANHANG AAANHANG 0.0000000
## ATRACTOR ATRACTOR 0.0000000
## AWERKT AWERKT 0.0000000
## ABROM ABROM 0.0000000
## ALEVEN ALEVEN 0.0000000
## APERSONG APERSONG 0.0000000
## AGEZONG AGEZONG 0.0000000
## AWAOREG AWAOREG 0.0000000
## AZEILPL AZEILPL 0.0000000
## APLEZIER APLEZIER 0.0000000
## AFIETS AFIETS 0.0000000
## AINBOED AINBOED 0.0000000
## ABYSTAND ABYSTAND 0.0000000Las variables PPERSAUT y MKOOPKLA son las más importantes
probs.test <- predict(bc, Caravan.test, n.trees = 1000, type = "response")
pred.test <- ifelse(probs.test > 0.2, 1, 0)
matriz <- confusionMatrix(as.factor(pred.test), as.factor(Caravan.test$Purchase))
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 4495 277
## 1 38 12
##
## Accuracy : 0.9347
## 95% CI : (0.9273, 0.9415)
## No Information Rate : 0.9401
## P-Value [Acc > NIR] : 0.9446
##
## Kappa : 0.0541
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.99162
## Specificity : 0.04152
## Pos Pred Value : 0.94195
## Neg Pred Value : 0.24000
## Prevalence : 0.94007
## Detection Rate : 0.93219
## Detection Prevalence : 0.98963
## Balanced Accuracy : 0.51657
##
## 'Positive' Class : 0
## Vemos que la predicción de personas que en realidad harán la compra tiene una precisión del 3.46%
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type
## == : prediction from a rank-deficient fit may be misleadingpred.test2 <- ifelse(probs.test > 0.2, 1, 0)
matriz <- confusionMatrix(as.factor(pred.test2), as.factor(Caravan.test$Purchase))
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 4495 277
## 1 38 12
##
## Accuracy : 0.9347
## 95% CI : (0.9273, 0.9415)
## No Information Rate : 0.9401
## P-Value [Acc > NIR] : 0.9446
##
## Kappa : 0.0541
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.99162
## Specificity : 0.04152
## Pos Pred Value : 0.94195
## Neg Pred Value : 0.24000
## Prevalence : 0.94007
## Detection Rate : 0.93219
## Detection Prevalence : 0.98963
## Balanced Accuracy : 0.51657
##
## 'Positive' Class : 0
## Vemos que usando una regresión logística la precisión del modelo a la hora de predecir personas que realmente realizaron la compra es de un 3.46%
12
Se usará la base de datos ‘Smarket’, el cual representa los porcentajes diarios de rendimiento para el índice bursátil S&P 500 entre 2001 y 2005. En este se intentará predecir la variable ‘Direction’
train <- sample(nrow(Weekly), nrow(Weekly)*0.7)
Weekly$Direction <- ifelse(Weekly$Direction == "Up", 1, 0)
Weekly.train <- Weekly[train, ]
Weekly.test <- Weekly[-train, ]- Aplicando boosting
bf <- gbm(Direction ~ . - Year - Today, data = Weekly.train, distribution = "bernoulli", n.trees = 5000)
bprobs <- predict(bf, newdata = Weekly.test, n.trees = 5000)
bpred <- ifelse(bprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(bpred), as.factor(Weekly.test$Direction))
matriz$table## Reference
## Prediction 0 1
## 0 84 100
## 1 54 89## [1] "Accuracy: 0.5291 Sensitivity: 0.6087 Specificity: 0.4709"- Aplicando bagging
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?bagprobs <- predict(bagf, newdata = Weekly.test)
bagpred <- ifelse(bagprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(bagpred), as.factor(Weekly.test$Direction))
matriz$table## Reference
## Prediction 0 1
## 0 56 57
## 1 82 132## [1] "Accuracy: 0.5749 Sensitivity: 0.4058 Specificity: 0.6984"- Aplicando random Forest
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?rfprobs <- predict(rff, newdata = Weekly.test)
rfpred <- ifelse(rfprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(rfpred), as.factor(Weekly.test$Direction))
matriz$table## Reference
## Prediction 0 1
## 0 58 54
## 1 80 135## [1] "Accuracy: 0.5902 Sensitivity: 0.4203 Specificity: 0.7143"- Regresión logística
logf <- glm(Direction ~ . - Year - Today, data = Weekly.train, family = "binomial")
logprobs <- predict(logf, newdata = Weekly.test, type = "response")
logpred <- ifelse(logprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(logpred), as.factor(Weekly.test$Direction))
matriz$table## Reference
## Prediction 0 1
## 0 7 8
## 1 131 181## [1] "Accuracy: 0.5749 Sensitivity: 0.0507 Specificity: 0.9577"- Regresión lineal
lmf <- lm(Direction ~ . - Year - Today, data = Weekly.train)
lmprob <- predict(lmf, Weekly.test)
lmpred <- ifelse(lmprob > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(lmpred), as.factor(Weekly.test$Direction))
matriz$table## Reference
## Prediction 0 1
## 0 7 8
## 1 131 181## [1] "Accuracy: 0.5749 Sensitivity: 0.0507 Specificity: 0.9577"En general no se encuentran grandes diferencias entre ellos, pero de estos propuestos, el que mejor resultados da es el bagging con un accuraccy de 54.74%
Capitulo 9
4
## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.gen <- function(x, y) {
x^2 + y^2 <= 1
}
df <- df %>%
rename(Var1 = X1, Var2 = X2) %>%
mutate(Clases = ifelse(gen(Var1, Var2),
'Clase 1',
'Clase 2'),
Clases = factor(Clases))
inTrain <- sample(nrow(df), nrow(df)*0.7)
train <- df[inTrain,]
test <- df[-inTrain,]
theme_set(theme_bw())
ggplot(df, aes(Var1, Var2, col = Clases)) +
geom_point(size = 2)
- Clasificador de soporte vectorial
svm1 <- svm(Clases ~ ., data = train,
kernel = 'linear',
scale = FALSE, cost = 10)
plot(svm1, data = train)
## Reference
## Prediction Clase 1 Clase 2
## Clase 1 8 10
## Clase 2 4 8## [1] "Accuracy: 0.5333 Sensitivity: 0.6667 Specificity: 0.4444"
- SVM con kernel polinomial
svm2 <- svm(Clases ~ ., data = train,
kernel = 'polynomial', degree = 2,
scale = FALSE, cost = 1)
plot(svm2, train)
## Reference
## Prediction Clase 1 Clase 2
## Clase 1 12 1
## Clase 2 0 17## [1] "Accuracy: 0.9667 Sensitivity: 1 Specificity: 0.9444"
- SVM con kernel radial
## Reference
## Prediction Clase 1 Clase 2
## Clase 1 11 2
## Clase 2 1 16## [1] "Accuracy: 0.9 Sensitivity: 0.9167 Specificity: 0.8889"
Podemos apreciar que, tal como lo decía el enunciado, tanto la SVM de kernel radial como polinomial tienen un mejor rendimiento que el modelo con kernel lineal. El modelo que tuvo un mejor rendimiento en este experimento fue la SVM polinomial de grado 2, con una precisión del 100% con el set de entrenamiento.
5
df <- data.frame(y,x1,x2)
df <- df %>%
mutate(y = factor(y ))
theme_set(theme_bw())
ggplot(df, aes(x=x1, y=x2, col=y)) +
geom_point()
##
## Call:
## glm(formula = y ~ ., family = "binomial", data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3452 -1.0996 -0.9886 1.1824 1.3738
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.038121 0.090234 -0.422 0.67269
## x1 0.008842 0.308455 0.029 0.97713
## x2 0.857039 0.309664 2.768 0.00565 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 692.86 on 499 degrees of freedom
## Residual deviance: 685.07 on 497 degrees of freedom
## AIC: 691.07
##
## Number of Fisher Scoring iterations: 4probs <- predict(logreg_fit, data = train, type = "response")
preds <- rep(0, 500)
preds[probs > 0.5] <- 1
train['preds'] <- predstheme_set(theme_bw())
train <- train %>%
mutate(preds = factor(preds))
ggplot(train, aes(x=x1, y=x2, col=preds)) +
geom_point()
##
## Call:
## glm(formula = y ~ poly(x1, 2) + x2, family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2464 -0.7807 -0.5253 0.7573 1.8684
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.05209 0.11049 0.471 0.637
## poly(x1, 2)1 0.16859 2.62143 0.064 0.949
## poly(x1, 2)2 32.04520 3.01261 10.637 <2e-16 ***
## x2 0.88434 0.36713 2.409 0.016 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 692.86 on 499 degrees of freedom
## Residual deviance: 516.82 on 496 degrees of freedom
## AIC: 524.82
##
## Number of Fisher Scoring iterations: 4probs <- predict(lm.fit, data = train, type = "response")
preds2 <- rep(0, 500)
preds2[probs > 0.5] <- 1theme_set(theme_bw())
train <- train %>%
mutate(preds2 = factor(preds2))
ggplot(train, aes(x=x1, y=x2, col=preds2)) +
geom_point()
theme_set(theme_bw())
train['preds3'] <- preds3
ggplot(train, aes(x=x1, y=x2, col=preds3)) +
geom_point()
svmr <- svm(y ~ x1 + x2, data = train,
kernel = 'radial',
scale = FALSE)
predsr <- predict(svmr, train)theme_set(theme_bw())
train['predsr'] <- predsr
ggplot(train, aes(x=x1, y=x2, col=predsr)) +
geom_point()
Podemos ver a través de estos experimentos que las mquinas de soporte vectorial con metodos no lineas en datos que no tienen “perimetro” lineal, son muy efectivas, también si tienen perímetro lineal las svm con método lineal también son muy efectivas en estos casos.
6
x1=runif (1000) -0.5
x2=runif (1000) -0.5
temp <- x1-x2
y <- rep(0, 1000)
for (i in 1:1000) {
if (temp[i] > 0.1) {
y[i] <- 1
}
else if (temp[i] < -0.1) {
y[i] <- 0
}
else {
y[i] <- runif(1) > 0.5
}
}
df2 <- data.frame(x1,x2,y)
df2$y <- as.factor(df2$y)
ggplot(df2, aes(x=x1,y=x2,col=y)) +
geom_point()
svm1 <- tune(svm, y ~ ., data = df2, kernel = "linear", ranges = list(cost = c(0.0001, 0.001, 0.01, 1, 5, 10, 100, 1000, 10000)))
summary(svm1)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 0.001
##
## - best performance: 0.086
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-04 0.480 0.03915780
## 2 1e-03 0.086 0.04427189
## 3 1e-02 0.092 0.03910101
## 4 1e+00 0.089 0.03984693
## 5 5e+00 0.092 0.04417138
## 6 1e+01 0.092 0.04237400
## 7 1e+02 0.089 0.03813718
## 8 1e+03 0.089 0.03813718
## 9 1e+04 0.091 0.04040077Podemos ver que un costo de 0.01 es el que nos dá el menor error
## cost misclass
## 1 1e-04 480
## 2 1e-03 86
## 3 1e-02 92
## 4 1e+00 89
## 5 5e+00 92
## 6 1e+01 92
## 7 1e+02 89
## 8 1e+03 89
## 9 1e+04 91Podemos ver que con un costo de 0.01 se clasificaron mal 85 observaciones.
set.seed(5)
x1_test=runif (1000) -0.5
x2_test=runif (1000) -0.5
temp_test <- x1_test-x2_test
y_test <- rep(0, 1000)
for (i in 1:1000) {
if (temp_test[i] > 0.1) {
y_test[i] <- 1
}
else if (temp_test[i] < -0.1) {
y_test[i] <- 0
}
else {
y_test[i] <- runif(1) > 0.5
}
}
df3 <- data.frame(x1_test,x2_test,y_test)
df3$y_test <- as.factor(df3$y_test)
ggplot(df3, aes(x=x1_test,y=x2_test,col=y_test)) +
geom_point()
names(df3) <- c("x1","x2","y")
costs <- c(0.0001, 0.001, 0.01, 1, 5, 10, 100, 1000, 10000)
test.err <- rep(NA, length(costs))
for (i in 1:length(costs)) {
svm.fit <- svm(y ~ ., data = df2, kernel = "linear", cost = costs[i])
pred <- predict(svm.fit, df3)
test.err[i] <- sum(pred != df3$y)
}
data.frame(cost = costs, misclass = test.err)## cost misclass
## 1 1e-04 486
## 2 1e-03 91
## 3 1e-02 92
## 4 1e+00 95
## 5 5e+00 94
## 6 1e+01 96
## 7 1e+02 96
## 8 1e+03 96
## 9 1e+04 96Podemos ver que en este caso el mejor costo fue dado por \(c=1\), con 92 errores a la hora de predicción
Podemos apreciar que aunque el costo en ambos casos fue diferente, realmente no difiere demasiado, ya que a partir de un costo de 0.001 no se presentan grandes cambios en los errores,
#### 7
Auto <- Auto[,-10] %>%
mutate(highmpg = as.integer(mpg > median(mpg))) %>%
mutate(highmpg = factor(highmpg))
head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name highmpg
## 1 chevrolet chevelle malibu 0
## 2 buick skylark 320 0
## 3 plymouth satellite 0
## 4 amc rebel sst 0
## 5 ford torino 0
## 6 ford galaxie 500 0svml <- tune(svm, highmpg ~ ., data = Auto, kernel = "linear", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000)))
summary(svml)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 1
##
## - best performance: 0.01275641
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.09205128 0.05715898
## 2 1e-02 0.07416667 0.04448352
## 3 1e-01 0.04865385 0.03720214
## 4 1e+00 0.01275641 0.01808165
## 5 5e+00 0.02044872 0.02020886
## 6 1e+01 0.02551282 0.01709615
## 7 1e+02 0.03826923 0.01345323
## 8 1e+03 0.03826923 0.01345323Se puede apreciar como con un costo de 1 se tiene el menor error.
- Kernel polinomial
svmk <- tune(svm, highmpg ~ ., data = Auto, kernel = "polynomial", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000), degree = c(2, 3, 4)))
summary(svmk)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost degree
## 1000 2
##
## - best performance: 0.2472436
##
## - Detailed performance results:
## cost degree error dispersion
## 1 1e-03 2 0.5635256 0.03971883
## 2 1e-02 2 0.5635256 0.03971883
## 3 1e-01 2 0.5635256 0.03971883
## 4 1e+00 2 0.5635256 0.03971883
## 5 5e+00 2 0.5635256 0.03971883
## 6 1e+01 2 0.5251282 0.11987686
## 7 1e+02 2 0.2959615 0.06069283
## 8 1e+03 2 0.2472436 0.05930020
## 9 1e-03 3 0.5635256 0.03971883
## 10 1e-02 3 0.5635256 0.03971883
## 11 1e-01 3 0.5635256 0.03971883
## 12 1e+00 3 0.5635256 0.03971883
## 13 5e+00 3 0.5635256 0.03971883
## 14 1e+01 3 0.5635256 0.03971883
## 15 1e+02 3 0.3341667 0.08484496
## 16 1e+03 3 0.2551282 0.06145828
## 17 1e-03 4 0.5635256 0.03971883
## 18 1e-02 4 0.5635256 0.03971883
## 19 1e-01 4 0.5635256 0.03971883
## 20 1e+00 4 0.5635256 0.03971883
## 21 5e+00 4 0.5635256 0.03971883
## 22 1e+01 4 0.5635256 0.03971883
## 23 1e+02 4 0.5635256 0.03971883
## 24 1e+03 4 0.5302564 0.09378263Aquí vemos que el menor error se da con un costo de 1000 de grado 2
- Kernel radial
svmr <- tune(svm, highmpg ~ ., data = Auto, kernel = "radial", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000), gamma = c(0.01, 0.1, 1, 5, 10, 100)))
summary(svmr)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost gamma
## 100 0.01
##
## - best performance: 0.01532051
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 1e-03 1e-02 0.59147436 0.06357909
## 2 1e-02 1e-02 0.59147436 0.06357909
## 3 1e-01 1e-02 0.08935897 0.03875505
## 4 1e+00 1e-02 0.07647436 0.03804385
## 5 5e+00 1e-02 0.05608974 0.03963690
## 6 1e+01 1e-02 0.03314103 0.02424635
## 7 1e+02 1e-02 0.01532051 0.01788871
## 8 1e+03 1e-02 0.03051282 0.02878864
## 9 1e-03 1e-01 0.59147436 0.06357909
## 10 1e-02 1e-01 0.24179487 0.10616790
## 11 1e-01 1e-01 0.08166667 0.03582639
## 12 1e+00 1e-01 0.05608974 0.03963690
## 13 5e+00 1e-01 0.02801282 0.02509397
## 14 1e+01 1e-01 0.02288462 0.02505026
## 15 1e+02 1e-01 0.03307692 0.03404377
## 16 1e+03 1e-01 0.03307692 0.03404377
## 17 1e-03 1e+00 0.59147436 0.06357909
## 18 1e-02 1e+00 0.59147436 0.06357909
## 19 1e-01 1e+00 0.59147436 0.06357909
## 20 1e+00 1e+00 0.06365385 0.04003438
## 21 5e+00 1e+00 0.06615385 0.04478910
## 22 1e+01 1e+00 0.06615385 0.04478910
## 23 1e+02 1e+00 0.06615385 0.04478910
## 24 1e+03 1e+00 0.06615385 0.04478910
## 25 1e-03 5e+00 0.59147436 0.06357909
## 26 1e-02 5e+00 0.59147436 0.06357909
## 27 1e-01 5e+00 0.59147436 0.06357909
## 28 1e+00 5e+00 0.52775641 0.07143278
## 29 5e+00 5e+00 0.52262821 0.07036637
## 30 1e+01 5e+00 0.52262821 0.07036637
## 31 1e+02 5e+00 0.52262821 0.07036637
## 32 1e+03 5e+00 0.52262821 0.07036637
## 33 1e-03 1e+01 0.59147436 0.06357909
## 34 1e-02 1e+01 0.59147436 0.06357909
## 35 1e-01 1e+01 0.59147436 0.06357909
## 36 1e+00 1e+01 0.55583333 0.05641062
## 37 5e+00 1e+01 0.54820513 0.06563816
## 38 1e+01 1e+01 0.54820513 0.06563816
## 39 1e+02 1e+01 0.54820513 0.06563816
## 40 1e+03 1e+01 0.54820513 0.06563816
## 41 1e-03 1e+02 0.59147436 0.06357909
## 42 1e-02 1e+02 0.59147436 0.06357909
## 43 1e-01 1e+02 0.59147436 0.06357909
## 44 1e+00 1e+02 0.59147436 0.06357909
## 45 5e+00 1e+02 0.59147436 0.06357909
## 46 1e+01 1e+02 0.59147436 0.06357909
## 47 1e+02 1e+02 0.59147436 0.06357909
## 48 1e+03 1e+02 0.59147436 0.06357909Con un costo de 100 y un gama de 0.01 se obtuvo el menor error en el modelo
svm.linear <- svm(highmpg ~ ., data = Auto, kernel = "linear", cost = 1)
svm.poly <- svm(highmpg ~ ., data = Auto, kernel = "polynomial", cost = 100, degree = 2)
svm.radial <- svm(highmpg ~ ., data = Auto, kernel = "radial", cost = 100, gamma = 0.01)- Lineal
- Polinomial
- Radial
8
data('OJ')
inTrain <- sample(nrow(OJ), 800, replace = FALSE)
training <- OJ[inTrain,]
testing <- OJ[-inTrain,]svm_linear <- svm(Purchase ~ ., data = training,
kernel = 'linear',
cost = 0.01)
summary(svm_linear)##
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "linear",
## cost = 0.01)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.01
##
## Number of Support Vectors: 445
##
## ( 222 223 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMPodemos apreciar que se tienen 440 vectores de soporte de los cuales 219 pertenecen a CH y 221 a MM
train.pred <- predict(svm_linear, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 411 76
## MM 64 249## [1] "Accuracy: 0.825 Sensitivity: 0.8653 Specificity: 0.7662"## error
## 0.175test.pred <- predict(svm_linear, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 158 27
## MM 20 65## [1] "Accuracy: 0.8259 Sensitivity: 0.8876 Specificity: 0.7065"## error
## 0.1740741Podemos ver que en ambos casos el error es del 17%
to = tune(svm, Purchase ~ ., data = training, kernel = "linear", ranges = list(cost = 10^seq(-2,
1, by = 0.25)))
summary(to)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 1.778279
##
## - best performance: 0.16625
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01000000 0.18000 0.04495368
## 2 0.01778279 0.17500 0.03864008
## 3 0.03162278 0.17250 0.04518481
## 4 0.05623413 0.17000 0.04297932
## 5 0.10000000 0.17250 0.04440971
## 6 0.17782794 0.17125 0.04566256
## 7 0.31622777 0.17625 0.04619178
## 8 0.56234133 0.17125 0.04450733
## 9 1.00000000 0.17000 0.04721405
## 10 1.77827941 0.16625 0.04678927
## 11 3.16227766 0.16750 0.04937104
## 12 5.62341325 0.16750 0.05075814
## 13 10.00000000 0.16875 0.05008673svm.linear <- svm(Purchase ~ ., kernel = "linear", data = training, cost = to$best.parameters$cost)
train.pred = predict(svm.linear, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 415 72
## MM 60 253## [1] "Accuracy: 0.835 Sensitivity: 0.8737 Specificity: 0.7785"## error
## 0.165test.pred <- predict(svm_linear, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 158 27
## MM 20 65## [1] "Accuracy: 0.8259 Sensitivity: 0.8876 Specificity: 0.7065"## error
## 0.1740741Podemos ver que el error de entrenamiento baja un poco, pero el de prueba se mantiene
##
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "radial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 372
##
## ( 190 182 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMtrain.pred = predict(svm.radial, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 428 78
## MM 47 247## [1] "Accuracy: 0.8438 Sensitivity: 0.9011 Specificity: 0.76"## error
## 0.15625test.pred = predict(svm.radial, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 159 34
## MM 19 58## [1] "Accuracy: 0.8037 Sensitivity: 0.8933 Specificity: 0.6304"## error
## 0.1962963to = tune(svm, Purchase ~ ., data = training, kernel = "radial", ranges = list(cost = 10^seq(-2,
1, by = 0.25)))
summary(to)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 3.162278
##
## - best performance: 0.17625
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01000000 0.40625 0.04903584
## 2 0.01778279 0.40625 0.04903584
## 3 0.03162278 0.32625 0.06440163
## 4 0.05623413 0.20250 0.05263871
## 5 0.10000000 0.18500 0.04281744
## 6 0.17782794 0.18375 0.04291869
## 7 0.31622777 0.17750 0.04322101
## 8 0.56234133 0.17750 0.05027701
## 9 1.00000000 0.17875 0.05369991
## 10 1.77827941 0.17750 0.05163978
## 11 3.16227766 0.17625 0.05350558
## 12 5.62341325 0.17625 0.05015601
## 13 10.00000000 0.17750 0.05130248svm.radial = svm(Purchase ~ ., data = training, kernel = "radial", cost = to$best.parameters$cost)
train.pred = predict(svm.radial, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 430 69
## MM 45 256## [1] "Accuracy: 0.8575 Sensitivity: 0.9053 Specificity: 0.7877"## error
## 0.1425test.pred = predict(svm.radial, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 160 35
## MM 18 57## [1] "Accuracy: 0.8037 Sensitivity: 0.8989 Specificity: 0.6196"## error
## 0.1962963##
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "poly",
## degree = 2)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 1
## degree: 2
## coef.0: 0
##
## Number of Support Vectors: 459
##
## ( 235 224 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMtrain.pred = predict(svm.poly, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 430 107
## MM 45 218## [1] "Accuracy: 0.81 Sensitivity: 0.9053 Specificity: 0.6708"## error
## 0.19test.pred = predict(svm.poly, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 171 46
## MM 7 46## [1] "Accuracy: 0.8037 Sensitivity: 0.9607 Specificity: 0.5"## error
## 0.1962963to = tune(svm, Purchase ~ ., data = training, kernel = "poly", degree = 2,
ranges = list(cost = 10^seq(-2, 1, by = 0.25)))
summary(to)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 10
##
## - best performance: 0.185
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01000000 0.40375 0.05894029
## 2 0.01778279 0.37750 0.06449591
## 3 0.03162278 0.36875 0.06594663
## 4 0.05623413 0.34625 0.06237487
## 5 0.10000000 0.31750 0.06487167
## 6 0.17782794 0.24125 0.06292908
## 7 0.31622777 0.21750 0.04609772
## 8 0.56234133 0.21125 0.04619178
## 9 1.00000000 0.20875 0.04641674
## 10 1.77827941 0.19375 0.04649149
## 11 3.16227766 0.19000 0.03809710
## 12 5.62341325 0.18750 0.04124790
## 13 10.00000000 0.18500 0.03855011svm.poly = svm(Purchase ~ ., data = training, kernel = "poly", degree = 2, cost = to$best.parameters$cost)
train.pred = predict(svm.poly, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 433 74
## MM 42 251## [1] "Accuracy: 0.855 Sensitivity: 0.9116 Specificity: 0.7723"## error
## 0.145test.pred = predict(svm.poly, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table## Reference
## Prediction CH MM
## CH 164 37
## MM 14 55## [1] "Accuracy: 0.8111 Sensitivity: 0.9213 Specificity: 0.5978"## error
## 0.1888889Podemos ver como el polinomial usando el mejor costo es el que ofrece un mejor resultado de entrenamiento pero el lineal es el que tiene menor error de prueba
Ensayo
Medellín es una ciudad la cual debido a su crecimiento poblacional y posición geográfica, presenta en general unos altos índices de contaminación, en particular, en algunas épocas del año esta crisis ambiental se agudiza. Más del 80% de las personas que viven en áreas urbanas que monitorean la contaminación del aire están expuestas a niveles de calidad del aire que exceden los límites de la Organización Mundial de la Salud (OMS, 2016). Esto hace necesario que se busquen mecanismos que ayuden a mejorar la problemática de la calidad del aire.
Por un lado tenemos que debido a la topografía de Medellín, y en general del valle de Aburrá, en épocas del año en donde las condiciones climáticas son propicias para que el aire no pueda calentarse lo suficiente, y en consecuencia no pueda ascender para que esta sea arrastrado fuera del Valle por los vientos alisios. Desde este punto de vista, es importante profundizar en la investigación de los modelos predictivos meteorológicos, los cuales, al poder darnos alertas mucho más precisas, podrían usarse para tomar medidas oportunas y no especulativas de cara a afrontar el problema.
Para llevar esto a cabo se puede tomar como referencia al trabajo realizado por Y. Zheng et al. en el año 2015, en donde pronosticaron la lectura de una estación de monitoreo de la calidad del aire de las siguientes 48 horas utilizando un modelo basado en datos (data-driven). Este sistema cuenta con un predictor de regresión lineal para modelar los factores locales de la calidad del aire. Para modelar los factores globales se usa una red neuronal, el cual es un predictor espacial. También cuenta con un agregador dinámico el cual combina las predicciones de los predictores teniendo en cuenta también los datos metereológicos y por último cuenta con un predictor para capturar los cambios repentinos en la calidad del aire.
Pero para esto tal y como subrayan en el año 2017 Kaur et al. existe una necesidad por mejorar la toma de información, ya que en definitiva, de esta dependen todos los modelos que se planteen. En este apartado es importante que Medellín cuente con un sistema de medición en tiempo real en donde no solo se monitoreen algunas áreas de la ciudad, si no donde podamos tener un entendimiento de que es lo que está ocurriendo realmente en las diferentes calles de Medellín, para esto podríamos usar diferentes dispositivos ayudándonos del internet de las cosas. Uno de estos ejemplos lo presenta Wooden en el año 2019, en donde nos muestra como Intel y Bosch desarrollaron un sistema al que denominaron The Intel®-based Bosch Air Quality Micro Climate Monitoring System (MCMS), el cual sin ocupar mucho espacio, podría usarse en diferentes zonas de la ciudad, para saber no solamente, si en general en Medellín hay mucha contaminación, si no en qué lugares específicamente es donde más contaminación se está presentando, en que calles en específico y a que horas del día, para poder tomar más medidas en estos lugares, ya que la contaminación que allí se genera termina teniendo un impacto en una zona mucho más amplia y afecta directamente a una gran parte de la población. Esto no solo ayudaría a tomar medidas, si no que permitiría evaluarlas para saber si realmente estas eran el problema principal.
- Referencias:
Organización Mundial de la Salud. (2016). Ambient Air Pollution: A global assessment of exposure and burden of disease. World Health Organization, 1–131. Retrieved from www.who.int.org
Y. Zheng, X. Yi, M. Li, R. Li, Z. Shan, E. Chang and T. Li, “Forecasting Fine-Grained Air Quality Based on Big Data,” Proceeding KDD ’15 Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 2267-2276, August 10, 2015.
Kaur, Gaganjot & Gao, Jerry & Chiao, Sen & Lu, Shengqiang & Xie, Gang. (2017). Air Quality Prediction: Big data and Machine Learning Approaches.
Wooden, A. (2019). Tackling air pollution with big data analytics. Intel. Recuperado de https://www.intel.co.uk